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SYSTEM  AVAILABILITY:   TIME  DEPENDENCE  AND  STATISTICAL  INFERENCE 

BY  (SEMI) NON-PARAMETRIC  METHODS. 

D.P.  Gaver 
P . A .  Jacobs 


INTRODUCTION. 

Component  or  system  avai labi 1 ity  generally  refers  to  (a)  the 
obability  that  the  item  of  concern  is  operable  or  "up"  and  mission- 
pable  at  a  point  or  during  a  period  of  time,  or  to  (b)  the  fraction 
time,  or  of  demands  if  the  item  is  on  standby,  for  which  the  system 
up.   Applied  to  large  systems  such  as  entire  U.S.  Navy  ships,  system 
ai lability  is  referred  to  as  read  i  ness .   At  a  different  level,  the 
ai lability  and  fault-tolerance  of  computer  systems,  both  government 

d  commercial,  is  of  concern,  and  has  recently  been  extensively 

•     •  1 
del  led  and  simulated;  see  Goyal  ,  Lavenberg  and  Trivedi  .   The 

ai lability  concept  is  also  relevant  when  discussing  the  safety  and 

oductivity  of  commercial  nuclear  power  plants;  in  that  arena  it  is 

antif ied  by  probabilistic  risk  assessment  (PRA) .   Related  finite 

2 
ate  stochastic  models  also  occur  in  medical  studies;  c.f .  Cox  , 

cobs  . 

Component  or  system  availability  is  influenced  both  by  the 

lerent  f ai 1 ure-proneness  of  the  item  and  by  the  time  and  resources  it 

kes  to  restore  a  failed  item  to  service.   Times  to  failure  or  "up 

nes"and  to  restoration  or  "down"  times  may  vary  considerably,  and  not 

:essarily  independently,  depending  upon  the  mode  of  failure,  the  time 

quired  to  diagnose  the  failure  including  the  access  to  (including 


compet it,  ion  for)  diagnostic  equipment  and  human  skills,  the 
availability  of  spare  parts  (logistic  delays),  and  other  factors.    It 
is  also  quite  conceivable  that  the  quality  of  the  repair  activity 
influences  future  times  to  failure.   This  effect  is  not  recognized  by 
the  usual  models;  but  see  Gaver   (pp.  775-800),  and  Thomas,  Jacobs  and 
Gaver  .    Incorporation  of  "availability  growth"  in  availability 
situations  is  a  practical  issue  that  is  not  frequently  modelled;  see 
our  Section  2  below,  however. 

The  present  paper  addresses  the  assessment  of  simple  system 
availability  when  there  is  concern  about  (a)  time-dependence ,  so 
demands  for  system  performance  are  not  necessarily  when  the  system  is 
in  "steady-state,"  as  is  often  assumed,  and  when  (b)  information  about 
system  failures  and  repairs  is  in  the  form  of  observed  data  so 
questions  of  statistical  influence  arise.   The  methods  and  models 
involved  lean  towards  the  semi -paramet r i c  or  non-parametric;  in 
particular  we  employ  the  empirical  Laplace  transform  in  the  time- 
dependent  scenarios  of  interest;  non-parametric  estimation  in  assumed 

steady-state  situations  has  been  studied  by  Gaver  and  Chu  .   Other 

7  S  Q 

investigations,  e.g.  Cox  ,  Gross  and  Harris   and  Ascher  and  Feingold 

have  been  overwhelmingly  concerned  with  estimation  in  presumably  we  1 1 - 

specified  stochastic  (queueing,  Markov)  models;  considerations  of  model 

mis-specification  are  seldom  broached.   Here  we  propose  analytically 

simple  approximations  to  time-dependent  system  behavior,  and  assess  the 

effects  of  model  specification  ("up"  and  "down"  time  dependence)  upon 

rates  of  approach  to  a  long-run  steady  state  as  the  latter  are 

estimated  from  available  data  (assumed  to  be  a  random  sample) . 

Uncertainty  assessments  (confidence  limits  and  standard  errors)  are 


furnished.   More  elaborate  procedures  involving  Bayes  or  empirical 
Bayes  setups  that  permit  "strength  borrowing"  (in  John  Tukey's  phrase) 
are  not  addressed  here,  but  are  agenda  items. 

Our  paper's  plan  is  as  "follows.   Section  2  describes  a  selected 
group  of  probability  models  for  simple  system  availability;  no 
comprehensiveness  is  claimed.   Solutions  are  given  in  terms  of  Laplace 
transforms,  all  of  which  are  rendered  immediately  i nterpretabl e  in 
terms  of  the  random-t ime-of -demand  or  observation  paradigm,  described 
originally  and  applied  to  transform  inversion  in  Gaver   .   Section  3 
describes  and  "fits"  the  simple  exponent ial -approach-to-steady-state 
model  used  for  representing  time-dependent  behavior;  see  Odoni  and 

1  Q  11 

Roth    and  Gaver  and  Jacobs   where  such  an  idea  was  used  to  represent 
time-dependent  queueing  behavior.   Section  4  introduces  issues  of 
estimating  time-dependent  availability  where  only  statistical  data  is 
at  hand.   Section  5  presents  a  variety  of  numerical  illustrations. 


2.    TIME-DEPENDENT  AVAILABILITY:   MODELS,  TRANSFORMS,  AND  THE  RANDOM 
OBSERVATION  (ROBS)  INTERPRETATION. 

Consider  a  system  that,  is  active  or  operative  ("up")  until  failur 
occurs,  after  which  a  "down"  period  occurs,  at  the  termination  of  whic 
a  new  up  period  begins,  and  so  on.   Suppose  the  system  is  to  be 
responsive  to  a  randomly  appearing  demand;  what  is  the  probability  tha 
it  will  respond  properly,  or  be  available,  at  the  time  of  the  demand? 
In  the  military  arena  the  device  might  be  a  surveillance  (radar,  sonar 
or  communications  system.    In  the  commercial  nuclear  power  context  it 
might  be  a  coolant  pump. 

There  are  many  versions  of  this  problem.   The  simplest  and  most 
traditional  imagines  that  the  system  operates  continually  between 
failures,  and  failure  is  instantly  detected  upon  occurrence.   Up  times 
may  be  modelled  as  iid  random  variables,  either  directly  or  indirectly 
If  the  system  is  redundant  a  new  up  time  begins  at  the  instant  an  old 
one  terminates,  with  repair  off-line.   Down  times  may  also  be  iid. 
However,  there  are  many  plausible  exceptions  to  these  scenarios,  a  few 
of  which  are  considered  next. 

Here  are  various  sample  model  formulations;  the  list  is  by  no 
means  complete.   In  all  of  these  the  system  is  up  at  t  =  0; 
modifications  are  straightforward. 

Mode  1  1  :   up  times  are  iid  {Uj}  ,  with  df  Fy(x);  down  times  are  iid 
{D|} ,  d.f .  FD(y)  .   A(t)  =  P{System  up  at  t | up  time  starting  at  t  =  0}  . 
Then,  by  a  simple  backward  renewal  argument,  A(t)  satisfies  the 
integral  equation 


A(t)  =  Fu(t)  +     A(t-v)Fc(dv) 

J0 


(2.1) 


'here  Fc(t)  =    Fu*FD(t),  the  convolution.  Laplace-transforming, 


A(s)  = 


oo 


e"stA(t)dt 


1  -  Fu(s) 


(2.2) 


If  the  item  is  observed  or  demanded  at  random  time  T  ~  exp(s)  then 


/.OO 


E[A(T)]  = 


-st 


""Awdt  =  i :  lucU 


(2.3) 


an 


d  since  E [T]  =  l/s  =  t,     the  probability  of  being  up  upon  demand 


1  s 


-l 


a0(r)  =  E[A(T)]  = 


1  -  Fy(r-A) 
1  -  FVCr"1) 


(2.4) 


E  [U]        E  [U] 


E[C]    E[U]  +  E[D] 


as   r  — <■  oo 


(2.5) 


by  Tauber ian/Abe 1 i an  results;  see  Feller  (1966) .   One  can  view  (2.4)  as 
the  availability  under  random  demand  or  observation;  ROBS  for  short. 
To  emphasize  the  fact  that  exponential  random  observation  is  involved 
we  utilize  EROBS .   The  last  formula,  written  here  as 


a0(oo) 


E[U] 


E[U]  +  E[D] 


(2.6) 


is  the  widely-applied  long-run  availability.    It  may  well  be  inappro- 
priate in  many  of  the  contexts  in  which  it  is  applied. 
Note ,  though ,  that 


a0(r)  = 


1  -  Fu(r-1) 


-1> 


FC(0 


(2.7, 


has  a  definite  interpretation  for  al 1  r  .    It  is  usually  easily  compute 
and  interpreted  whenever  transforms  of  U  and  D  (or  C)  exist.   As  will 
be  seen,  it  can  be  estimated  from  data  under  many  circumstances — even 
when  the  model  selected  has  either  E [U]  or  E [D]   infinite  and  (2.6) 
becomes  uninteresting. 

Here  are  several  additional  models  whose  transforms  can  be 
directly  written  down,  and  directly  interpreted  under  EROBS  demands. 

Model  2:   The  initial  up  time  U0  has  df  Fy  ;  {Dj  +  Uj  =  Cj  , 
i  =  1,2,...}  are  i id  and  independent  of  U0 .   This  attempts  to  model  a 
situation  in  which  the  quality  of  maintenance  during  a  down  time 
affects  the  distribution  of  the  next  up  time,  so  these  times  are 
dependent;  the  influence  stops  at  that  point,  however.   Let  Ap(t)  = 
P{system  not  up  (down)  at  t | down  time  starts  at  t  =  0}.   Then 


AD(t)  =  FD(t)  +     AD(t-r)Fc(dr) 

J0 


and 


A(t)  = 


so 


1   if   U0  >  t, 


AD(t-U0)  =  1  -  AD(t-U0)   if   U0  <  t 


A(t)  =  FUo(*)  +    C1  "  AD(t-r)]Fu  (dr)  =  1 


AD(t-r)FU()(dr) 


Transforming  gives 


A(s) 


1  -  Fc(s)  -  FUq(s)  +  FUq(s)Fd(s) 
1  -  Fc(s) 


(2.8) 


from  which  the  observational  probability  aQ(r)  can  be  immediately 
written  down:   a0(7-)  =  r"  A(r"  ).   Models  that  propose  that  the  up  time 
can  positively  influence  the  duration  of  the  following  down  time  have 
been  formulated;  see  Gaver 

Model  3:    (Random  changepoint  model;  "reliability  growth").   Suppose 
the  system  f  ai  1  s  accord  ing  to  up  t  imes  {Ujl5  i  =  1,2,  ...I},  iid  Fj  (u  ) 
and  thereafter  {Uj2  ,  i  =  1  +  1,1+2,...},  iid  F2(u);  I  ~  Pj  is  the  change- 
point  and  is  associated  with  d iagnost i c- repai r  activity.   Let  FD(u) 
describe  the  iid  down  times  {Dj  ,  i  =  1,2,...}.   Suppose  the  system 
starts  operation  as  a  Type  1  up  time,  transitioning  after  a  random 
number  of  down  times  to  Type  2.   Let  Aj:(t)  denote  the  probability  that 
the  system  is  up  and  in  state  j  at  time  t  given  it  is  up  at  time  0, 
i.e.  t  is  contained  in  an  up  time  of  Type  j.   Then,  using  backward 
renewal  arguments , 


AnOO  =  F1LJ(t)  +  (1-Pl)     A11(t-r)Flc(dr) 

J0 


A22(t)  =    F2U(t)  + 


ft 


A22(t-r)F2c(dr)   , 


and 


t 

A12(t)       A22(t-r)G(dr)  , 
0 


(2.9) 


where 


oo 


G(r)  =  £  PrTcCO   ^ 
i=l 


as  before  FjC(r)  refers  to  the  df  of  a  cycle  time  Ujj  +  Dj  .   Transforminj 
immediately  produces  a  solution: 


An(s) 


1     1  -  F1LJ(s) 

s  1  -  (l-Pl)Flc(s) 


A22(s)  -  s 


1  1  -  F2U(s) 
-  F2C(s) 


(2.1( 


A12(s)  =  A22G(s)  , 


and 


O)  =  £  p^icC^))'  =  p(f,ic(s))  . 


this  last  being  the  generating  function  of  the  number  of  repairs  to 
transition  rv , I .   Assemble  to  obtain  the  transform  of  the  probability 
of  being  up  at  t: 


-st 


Ai(s)  =     e-"[An(t)  +  A12(t)]dt  =  An(s)  +  A12(s) 
J0 


1  ~  F1LJ(s) 
1  -  (l-Pl)Flc(s) 


+  p(Fic(s))  • 


1  -  F2U(s) 


2U 
F2C(S) 


(2.11 


Consequently    the    observational     probability 


a0(r)     =    E[A(T)]     = 


1     -     d-Pi)Flc(l/r)  ^    1C'  n    1    -    F 


F2U(1/Q 


F2C(l/r) 


(2.11 


E[U2] 


E[U2]     +    E[D] 


as 


oo 
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but  the  latter  limiting  formula  provides  no  information  about  the 
availability   at  early  times,  possibly  before  change  has  taken  place. 
The  EROBS  random  demand  formulation  at  least  allows  some  information  to 
be  obtained  in  a  very  direct  and  meaningful  manner  by  just  evaluating 
the  transform  itself  numerically. 

Model  4:    (Markovian  changepo i nts) .    It  is  possible  that  a  system 
alternates  slowly  between  failure  states,  either  temporarily  or 
forever.   The  random  motion  may  be  in  response  to  occasional  changes  in 
maintenance  practices,  to  debugging  (reliability  growth) ,  or  to  ageing 
(reliability  decay)  and  subsequent  replacement.   To  illustrate,  let 
there  be  just  two  failure  distribution  states,  as  in  Model  3,  but  let 
{pi!  ;  i  ,  j  =    1?2}-  be  the  transition  probabilities  of  a  Markov  chain  that 
governs  jumps  between  them.   Suppose  the  system  starts  in  state  1,  at 
the  beginning  of  an  up  time  U(  ~  Fjy(r),  i  =  1,2.   Then,  letting  Ajj(t) 
be  the  probability  the  system  is  up  and  in  state  j  at  time  t  given  it 
just  started  an  up  time  at  time  0  and  was  in  state  i , 


An(t)  =  F1LJ(t)  +  p 


n 


A11(t-r)Flc(dr)  +  p12     A21  (t-r)  F1C  (dr)  , 
0  J0 


•t  rt 

Ai2(t)  P12  A22(t-r)Flc(dr)     +    Pll  A12  (t-  r )  F1C  (dr)  , 

0  J0 


(2.13) 


A22(t)    =    F2U(t)    +    p22  A22(t-r)F2C(dr)    +    p21 

J0 


A12(t-r)F2C(dr)  , 


t 

A2i(t)  P21  An(t-r)F2C(dr)    +    p22 

0 


A21(t-r)F2C(dr) 


Transforming    yields    these    equations, 


An(s)     = slu(S)    +    PiiAu(s)Flc(s)     +    Pl2A21(s)Flc(s) 


A12(s)     =    Pl2A22(s)Flc(s)     +    PllA12(s)Flc(s) 


A22(s)    =    s2LjC         +    p22A22(s)F2C(s)    +    p21A12(s)F2C(s) 


A21(s)     =    p21A11(s)F2C(s)     +    p22A21(s)F2C(s)  . 


Solve     (2. 14, a)     and     (2.14,d)     simultaneously    "to    get 


Au(s) 


1    "    Fiu(s) 


1     -     Pll^lcO)      -     Pl2P2lF,lc(s)F,2c(s)/(1-P22^2c(s)) 


From    (2. 14b, c) , 


Ai2(s)     -     s 


1    -    F2U(s) 


l-p22F2C(s)-p21p12Flc(s)F2C(s)/^l-PllFlc(s)) 

/    P12Flc(s)     x 
Vl-PllFlc(s)^ 


It    follows    that     if    t    —    1/s, 


ai(r)     =    sAlt    =    s[An(s)     +    A12(s)] 


1  Fiu(s) 


[l-p11Flc(s)-p12p21Flc(s)F2C(s)/(l-p22F2C(s)J 

1    -    F2U(s)  


I 


l-p22F2C(s)-p12p21Flc(s)F2C(s)/(l-PllFlc(s)J 

(    Pl2^ic(s)     ) 
Vl-PllFlc(s)' 

10 


The  above  formula  is  manifestly  opaque  as  it  stands.   Tauberian/ 
ibel ian  results  present  the  long-run  behavior: 


all(oo)  =  lim  sAll(s)  =  ELCl]  +  ^/p,,)  K  [C2]  <2-18) 


md 


ai2(0O)    =    E[C2]     +    (S/p12)E[C1]  (2-19> 


so    the    sum 


/     n  /      n  /      x  E[U,]  E[U2] 

iCoo)     =    an(oo)     +    a12(oo)     =  A  + 


ECCi]     +    (P12/P21)E[C2]     T    E[C2]     +    (P21/P12)E[C1] 

(2.20) 


•epresents  the  long-run  probability  that  the  system  is  up.   Of  course 
1^(7-)  for  r    finite  is  available  from  (2.17)  and  potentially  provides 
:onsiderably  more  information. 

Many  other  such  models  can  be  constructed  and  "solved"  by 
transforms.    In  what  follows  we  illustrate  the  way  in  which  such 
transform  solutions  can  be  exploited  to  yield  time-dependent  and 
nferential  information.   The  intention  is  to  provide  simple  and 
?lexible  approximate  information  rather  than  to  utilize  more  exact,  but 
lIso  more  computationally  demanding,  procedures. 
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3.    TIME-DEPENDENT  AVAILABILITY  AT  FIXED  OBSERVATION  TIMES  (FOBS): 
SIMPLE  EXPONENTIAL  REPRESENTATIONS. 

It  is  familiar  from  the  theory  of  finite  Markov  chains  and  also 

from  renewal  theory  that  if  a  long-run  or  steady-state  condition  is 

reached  by  a  stochastic  model,  then  the  nature  of  the  approach  is  oftei 

essentially  exponential.   In  particular,  if,  in  Model  1,  U  ~  exp(A), 

D  ~  exp(>)  then  for  A(0)  =  1, 


A(t)  =  X-*-  +  ,-£-  e-(A+/i)\  (3.1) 


Such  a  time  dependence  is  exhibited  approximately  by  the  M/G/l  queue; 
see  Odoni  and  Roth   ;  Morse       ,  and  later  Kielson   ,  and  others  have 
assessed  the  relaxation  times  of  various  queueing  systems.   See  also 
Gaver  and  Jacobs    for  a  brief  discussion  of  the  time-dependent  M/G/l 
queue  utilizing  such  a  representation. 

It  is  we  1 1 -understood  that  natural  generalizations  of  the  Markov 
formulation  producing  (3.1)  can,  for  certain  models,  lead  to  damped 
oscillatory  approach  to  a  steady-state  value  a(oo)  ,  or  to  ultimate 
approach  that  is  exponential  modified  by  a  negative  power  of  t. 
However,  in  what  follows  we  shall  assume  that  the  "error"  at  time  t  is 
approximately  exponential: 


e(t)  =  A(t)  -  A(oo)  ~  ae_/?\  (3. 2, a) 


or,     more    elaborately, 


,  .        N     -(r+i#)t  ,  .        x     -(r-i0)t  /0    _    ,, 

~     (a1+ia2)ev  +     (aria2)ev  (3.2,b; 
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and  investigate  ways  of  assessing  parameter  values  from  (a)  perfectly 
specified  models,  wherein  the  distributions  of  U,  D  are  presumed  known, 
and  (b)  from  data.   This  section  deals  with  problem  (a). 

The  motivation  for  considering  simple  exponential  model  (3.2)  is 
the  desire  for  an  easily  comprehended  and  computed  assessment  of  time- 
dependent  availability  at  a  fixed  demand  or  observation  time 
(abbreviated  FOBS) .    In  many  cases  encountered  such  an  assessment  (or 
assessments)  of  a    and  /?  or  of  al  ,  a2  ,  r  and  9    can  provide  a  useful  sense 
of  the  behavior  of  A(t)  as  time  progresses  without  the  necessity  of  an 
extremely  time-consuming  and  computer- i ntens ive  transform  inversion 
procedure  or  of  a  symbolic  inversion  in  terms  of  polynomial  roots;  the 
latter  is  "explicit,"  but  hardly  comprehensible  in  general. 
Simulations  and  numerical  solutions  of  governing  integral  equations  are 
also  useful  approaches,  but  they  are  generally  more  computationally 
intensive  than  our  proposals  here. 

We  suggest  several  ways  of  matching  an  exponential  to  A(t)  by 
utilizing  the  transform  A(s) ,  presumed  given. 

3.1   Method  1:   Least  Squares,  Unweighted  and  Weighted  ("Tuned") 
Begin  by  considering  (3 . 2 , a) .    If 


e(t)    =    A(t)     -    A(oo)  ; 


(3.3) 


we    wi 


ish    to    represent    this    "error"     by    e(t)     =    ore 


-0t 


Cons  ider 


A(a,/?) 


f°°  \  f°°/  /?t\ 

(e(t)-e(t)  )2dt    =  e(t)-«e"P    )2dt     , 


(3.4) 
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the  integrated  squared  error;  the  object  is  to  minimize  A  by  choice  o- 
a    and  j3 .       Differentiation  on  a    easily  gives  the  optimizing  condition 


af 


=  2/? 


oo 


0 


e(t)e"/?tdt  =  2/?e(/?)  =  2/3 


A(/?)  - 


A(oo)" 
0 


(3.5) 


i.e.  aQ     is  evaluated  in  terms  of  the  known  Laplace  transform  of  A(t) 
and  of  e(t)  ,  evaluated  at  /? .   To  find  a  minimizing  f3     it  is  possible  tc 
proceed  by  squaring  the  integrand  of  (3.4) ,  discarding  the  e  (t)  part, 
and  substituting  aQ    for  a.   Simplification  leads  to 


■A(a0(/?),/?)  =  A(/?)  =  -2/?[e(/?)J 


(3.6) 


which  can  be  readily  searched  for  a  global  maximum,  (3Q  .    It  is  easily 
verified  analytically  that  if  U  ~  exp(A),  D  ~  exp(/i)  in  Model  1  then 
/?0  =  A  +  /z ,  as  it  should.    In  general  the  search  of  (3.6)  must  be 
conducted  numerically. 

The  straightforward  least-squares  procedure  can  be  tuned  towards 
different  time  ranges  by  appropriate  weighting.   We  may  seek  to 
minim  ize 


oo 


A(a,/?;w)  =     (e(t)-ae'^)2w(t;r)dt 
J0  v  ' 


(3.7) 


where 


w(t;r) 


1   near  r 

0   far  from  r 


(3.8) 
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len  the  minimizing  a    and  f3    values,  now  aQ(r)  ,  /?o(r)  j  should  lead  "to  an 
<ponential  approximation  that  performs  especially  well  near  r — 
Lthough  perhaps  less  well  further  away.   Again  it  is  convenient  to 
sleet  an  analytical  form  for  w(t;r)  that  is  compatible  with  the 
iplace  transform,   i.e.   involves  linear  combinations  of  exponentials. 
ri  example  is  the  density  function  of  exponential  order  statistics; 
bis  device  was  used  in  a  numerical  transform  inversion  scheme 
javer   ;  Gaver  and  Jacobs    mention  the  idea  in  the  present  context) . 
Here  is  a  specific  example  of  weighting.   Let 


w(t;r)  =  ae^VVl-e^)  =  3{^\-^\)  (3.9) 


here    a    =    -Jr- .       Then    minimization     leads    first    to 

n  or 


<9A 
da 


0 


(e(t)-ae-^)e-^.3(e-2/XV-e-3/iV)dt    =    0 


hich    gives 


V+W     -    e(/H3/0     -    a{^L-    -    ^)    =    0  (3.10) 


_  e(/?+2//)     -    e(/?+3/x) 

<*oW     -    i/(2/?+2/i)     -    1/(20+3/0 

r~^o    >>          -fajL.*    M     (2/^+2/x)  (2/?+3/z) 
=     [e(/?+2/0     -    e (0+3/0  J      ji ■  k^-11; 
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Squaring  and  omitting  "the  term  in  e  (t)  delivers 


-A(a0(/?),/?)  =  A(/?)  =  (2/?+2,i)(2/H3A0[e  (0+2/0 -e  (0+3/0  ] 


(3.12 


:edure  is  tuned  towards  r  =  -J^-     in  this  case.   A  sharper  tuning  can 

6/i 


which  can  be  easily  searched  for  a  global  maximum.   The  estimation 

proc 

be  accomplished  by  a  weight  function  that  concentrates  more  tightly 

around  r  than  does  (3.9);  the  latter  is  recognizably  the  density  of  th 

median  of  a  sample  of  size  3  from  an  exponential  df .    If  w  is  the 

density  of  the  mth  order  statistic  of  a  sample  of  n  exponentials  then 


its  mean  is  7= 


n  +  ^T  + 
n    n- 1 


+ 


n-m+1 


i  ln(l  -  ™-)  and  its  vari 


/< 


lance  u 


~  -h>        ,   m —  ..  N   if  m  is  proportional  to  n.   As   n  increases  the  density 
f    n(n-m+l) 

function  of  w  =  wn  approaches  a  Gaussian/Normal  density  with  variance 
decreasing  like  1/n ,  which  is  essentially  a  delta  function  at  r  as  n  - 

00  . 


3.2   Method  2:   Derivative  Matching:   Exploiting  Random  Observation 


Times  (ROBS) . 


-/?t 


Since  it  is  proposed  to  approximate  e(t)  =  A(t)-A(oo)  by  ae  ,  an 
since  the  transforms  of  both  sides  are  known  (or  can  be  estimated  fron 
data)  one  may  contemplate  determining  a    and  j3    as  the  solution  of 


e(r)  = 


00 


e"U/T[A(u)-A(oo)]  ^  =  ?A(l/r)-A(oo) 


00 


(I 


-u/r  -(3u    du        1 
e    ae    V  =  <*    IT^  5 


(3.13; 


that  is,  the  system  is  observed  at  a  random  exponential  time  having 
One  choice  for  r  is  t.   Obviously  any  a,     /?  so  found  must 


mean  r 
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mplicitly  depend  upon  r;  in  many  situations  for  which  the  exponential 
s  nearly  correct  the  dependency  of  /?  and  a    on  r  will  be  gentle.   Here 
re  two  derivative-like  prescriptions  for  finding  a(r),  /?(r): 
a)   Solve 


0 


I  +  g(r  +  A)  =  r(r  +  A)  =  l/e(r  +  A) 


(3.14) 


0 


J  +  g(r  -  A)  =  r(r  -  A)  =  l/e(r  -  A) 


or  a  =  a(r;A),  /?  =  /3(r;A). 

b)    Invert  and  analytically  differentiate  at  an  interesting  r-value 


r(r)  = 


A(l/r)  -  rA(oo) 


=  I  +  =-  r 
a  ~  a  ' 


/5 


r'(0  =  g  . 


/? 


(3.15) 


(3.16) 


Now  solve  simultaneously  for  /?(r)  ,  a(r)  .   A  good  diagnostic  step  is 
o  plot  q;(t)  and  /?(r)  vs.  r;  if  the  plots  are  nearly  constant  with  r 
hen  the  exponential  provides  a  useful  form;  if  the  change  is  gradual 
t  can  potentially  be  represented  by  an  empirical  function,  e.g.  a  low- 
rder  polynomial. 

It  is  also  possible  to  match  the  exponential  approximation  by  r- 
uning,  as  was  done  by  the  weighting  procedure  of  (3.7) .   Simply 
bserve  the  system  with  probability  density 


w(t;r)  =    c(n,m)(l-e"//t)m-1e"/iV(e"/it)n-r 


(3.17) 


t     is    convenient    to    take    n    odd    and    m    =     [n/2] 
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3.3   Method  3:   Matching  Transform  Means. 

The  following  procedure  represents  another  possibility  for 
assessing  /?  .   Put  a    =  1  -  A(oo)  and  write 

A(t)  =  A(oo)  +  (1  -  A(oo))P(t) 

and,     denoting    the    transform    of    P    by    p, 


r-OO 


0 


.-st 


se~slA(t)dt     =    sA(s)     =    A(oo)     +     (1     -    A(oo))sp(s) 


(3.19) 


Note    that     if 


A(t)     -    A(oo)  -0t 

1     -    A(oo) 


then  transformation  yields 


sp(s) 


1  +  /?(s) 


or 


so 


[i  +  /?(i)]sP(s)  =  i, 


0 


1  -  sp(s) 
P(s) 
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To    explicate    the    behavior    of    /?    for    large    t,     expand     (3.22)     in    powers    of 
s    for    s    smal 1 : 


P    = 


s[l     -    sA(s)] 
sA(s)     -    A(oo) 


sE[C]{Fc(s)     -    Fy(s)} 


E[C][1    -    Fu(s)]     -    E[U][1    -    Fc(s)] 


sE[C]{-s(E[C]-E[U])    +    %(E[C2]-E[U2])+o(s2)} 
E[C]{sE[U]     -    ^E[U2]}    -    E[U]{sE[C]     -    ^E  [C2]  }+o  (s2) 


Thus, 


(3.23) 


lim     /?(s)     = 
s— >-0 


2E[C]{E[C]     -    E[U]> 
E[C2]E[U]     -    E[C]E[U2] 


(3.24) 


It  is  anticipated  that  this  expression  is  useful  when  t  becomes  large, 
provided  a  simple  exponential  is  an  appropriate  approximation. 

We  will  evaluate  (3.22)  for  s  =  1/t,  t  being  the  time  of  interest 
It  is  also  possible  to  tune  this  procedure  to  a  time  t  as  before.    In 
particular,  setting 


oo 


0 


,    N  rA(t)  -  A(oo)i       I*01    ,.       s 
<*»0  {  1  -  A(oo)   }  dt  -    JQ   W(t;T)e 


■Pt 


(3.25) 


and  solving  for  /?  where,  for  example,  w  is  as  in  (3.9).   This  tactic 
will  be  employed  for  inferential  purposes  in  Section  5  with  w(t;r)  of 
the  form  in  (3.9) . 
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4.    TIME-DEPENDENT  APPROACH  CHARACTERISTICS  ASSESSED  FROM  DATA. 

In  this  section  we  examine  situations  reasonably  represented  as 
indicated  previously,  but  for  which  limited  information  on  the  up  and 
down  times  is  available:   only  samples  of  finite  sizes  nu  for  the  up 
times  and  nd  for  the  down  times  are  at  hand.   Furthermore  it  is  desirec 
to  avoid  using  particular  conventional  analytical  forms  for  the 
underlying  distributions;  instead  "non-parametric"  or  "distribution- 
free"  methods  will  be  invoked.   Finally,  the  availability  at  a  finite 
time  after  an  initial  moment  is  of  interest. 

4 . 1   Predicting  Availability  Assuming  EROBS :   Point  and  Interval 

Est  i mates . 

Suppose  we  wish  to  estimate  the  availability  of  a  Model  1-type 
system  at  a  random  (exponential)  time  with  mean  r,  i.e.  assuming  EROBS 
Data  are  available  on  the  up  times  and  down  times: 

ul  »  u2  '   •  •  •  >   unu  s 

dx  ,  d2  ,   .  .  .  ,  dn<j 
are  the  respective  observations;  although  the  order  of  observation  is 
Uj  ,  dj  ,  u2  ,  d2  ,   ...  etc.  ,  it  will  first  be  assumed  that  the  Uj  '  s  and  d:'i 
are  independently  sampled  from  fixed  distributions  Fy  and  FD ,  but  the 
latter  process  is  otherwise  unspecified. 

A  natural  non-parametric  procedure  is  to  define  the  sample  or 
empirical  transforms 


1   nu    -su;  -,   nd    -sd: 

*u(*0  =  ifc  £  «   '  .    *d(s)  =  rn  £  e   '  (4.1) 

i=l  d  i=l 


and  to  utilize  these  in  place  of  the  population  or  true  transforms  thai 
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appear  in  (2.4)  as  a  solution  to  the  problem  posed.   Thus  a  non- 
parametric  estimate  of  the  desired  availability  is 

ft  [A(T)]  =  a0(r)  =    J_^^>     .  (4.2) 

Notice  that  the  up  and  down  times  are  conveniently  assumed  independent 
here.   Procedures  to  deal  with  more  complex  models  can  and  should  be 
devised;  accommodation  to  possible  up  and  down  data  dependencies  is  a 
natural  step,  first  by  graphics  and  subsequently  by  model  fitting  and 
test  ing . 

There  are  various  options  for  estimating  the  sampling  variability 
of  the  estimate  (4.2)  given  the  basic  model.   A  classical  procedure, 
adopted  here,  is  to  first  note  that  under  our  assumptions  both  <^y(s) 
and  </>D(s)  are  asymptotically  Gaussian/Normal  by  the  central  limit 
theorem,  for  they  are  seen  to  be  (modelled  as)  averages  of  independent 
and  even  bounded  random  variables  e"s  ,  e"s  ;  s  =  1/r  >  0.   We  find 
easily  that  the  corresponding  random  variables  <^y(s)  ,  <^>q(s)  are 
approximately  distributed  as  follows: 

a    <    \  Jvm   Fu(2s)  -  (Fu(s))2A 

<t>u(s)  ~  n^Fu(s),  -y nn J  5 

(4.3) 

A>     (^           n(f  rO   FD(2s)  -  (FD(s))2x 
<pD(s)     ~  N(FD(s),  ^ J  . 

Since  a0(r)  is  a  probability,  confined  between  zero  and  one,  it  is 
reasonable  to  carry  out  further  asymptotics  on  a  transformation;  we 
choose  to  study  the  logistic  transformation 
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(r)  =    In 


a-oQ) 
l-a0(r) 


=  In 


1-*u(b) 
*u(s)(l-*D(s)) 


(4 


=  ln[l-^y(s)]  -  In  ^(s)  -  ln[l-*D(s)] 


Expand  in  Taylor  series  ("the  delta  method)  to  obtain 


E[£(r)]     ~     In 


In 


aoO) 
l-a0(r) 

a-p(r) 
l-a0(r) 


+ 


+ 


1-2F 


1    Var[y 


2  (a-^Fu)2     L'uJ     2  (i-fd)2 


l)(       1-2<^u       \„2 


sn/nd 


l(_i^u      U/nu  +   JM^ 


2r   ' 


(4.5} 


where  we  have  replaced  Fy  and  FD  by  their  estimates  <^>y  and  <pD     in  the 
correction  term.   Furthermore, 


Var[£(r)]  ~ 


(FuCi-Fu))2 

su/"u 

(^uCWu))*   (Wd) 


Var[^]  +     *   2  Var[y 
K1~t  d) 


+ 


sd/nd 


n2  ~  s£   ' 


(4.6 


We  have  put 


'"  =  n-7^T  £  (e_  '  -  *") 


and 


(4 


-,      nd  t    -sd:       v. 


Finally,  approximate  a-100%  confidence  limits  are  given  by 
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+  zft/2  •  ^  <  ln[a0(r)/(l-a0(r))]   <  I*    +  z^^s, 


(4.8) 


ind 


I  _  lnf Wu(lA) x  _  l(/    1-2^    x  4        1     4\  (4    Q) 

■      ~    inUu(l/r)(l-^D(l/r))J    ^((l-^)^)^  n"  +  (l-<f>D)2     "dj  '     ^     ^ 


vhile  s«  =  \sf    from  (4.6)  .   The  expression  (4.8)  is  "then  inverted  to 
provide  two-sided  confidence  limits  for  the  actual  availability.   A 
lumerical  example  is  provided  in  a  later  section. 

1.2   Predicting  Availability  Semi -Parametr i cal ly  at  Fixed  Time 
(According  to  FOBS) . 

Utilize  the  same  data  as  that  in  section  4.1  but  suppose  we  wish 

bo  estimate 


A(t)  ~  A(oo)  +  ae 


/?t 


(4.10) 


\fter  some  experimentation  it  has  been  found  that  Method  3  above  can  be 
nost  easily  and  effectively  adopted  to  estimate  0    and  a:        in  (3.18)  put 


A(oo)  = 


u  +  d 


(4.11) 


j-nd  estimate  (3    by  evaluating  the  empirical  version  of  (3.22).    In  the 
simplest  case  (Model  1),  the  empirical  version  of  p(s)  is  given  by 


sp(s)  = 


1-</>c(s)     u  +  d 


(4.12) 


u  +  d 
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so    "the    estimator    of     (3.22)     is 


3e(s)     = 


1     -    sp(s) 
POO 


(4.13! 


Sett  i  ng 


w(t;r)     =    3fi(e  -    e         ) 


in     (3.25)     results    in    the    equation 


3/i[p(2/0     -    p(3//)]     =    3/i 


(2/i    +    /?)(3//    +    /?) 


or 


/?2    +    5/i/?    +    6/i' 


p(2/i)     -    p(3j0 


=    0 


(4.14) 


a  quadratic  in  /?  .   The  positive  solution  to  (4.14),  /?q  ,  gives  another 
estimate  of  /? : 


£q    =    J  [-50    +    J/i2    +    4/x[p(2/i)     -    p^)]"1] 


(4.15) 


5  ~ 

where  /j  =  ^— .    In  practice  /?q  may  occasionally  be  negative  or  imagi- 
nary, in  which  case  a  reasonable  alternative  is  required.    In  our  simu 
lation,  tests  of  the  procedure  defaulted  to  an  estimate  appropriate  fo 
U,D  independent  and  exponentially  distributed,  i.e.  /?  =  (i  +  -Lj.   How- 
ever in  general  a  non-positive  J3    may  indicate  i  nappropr  iateness  of  the 
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exponent i al -approach  model  (4.10)  so  a  problem-dependent-  alternative 
must  be  sought . 

The  asymptotic  behavior  of  the  proposed  estimates  is  not  explored 
here.    In  the  following  section  we  describe  bootstrapping  (sample  re- 
use) procedures  and  results  for  accessing  sampling  variability  of 
estimates  of  A(t)  and  j3 . 
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5.   NUMERICAL  ILLUSTRATIONS  BY  SIMULATION. 

The  procedures  described  earlier  will  now  be  illustrated,  and  to 
some  degree  tested,  using  simulated  data. 

5.1   Predicting  Availability  under  EROBS :    Illustrations 

Refer  to  Section  (4.1)  for  the  basic  approach  to  be  utilized. 
Suppose  that  nu  =  nd  observations  are  available  on  i id  exponentially 
distributed  up  and  down  time  random  variables;  E [U]  =  1,  E [D]  =  4. 
These  data  are  to  be  processed  utilizing  expressions  (4.8)  and  (4.9) . 
To  do  so,  a  total  of  ns  —    500  independent  simulations  were  carried  out 
for  the  cases  nu  =  nd  =  n  =  10,  and  nu  =  nd  =:  n  =  25.   In  Table  5.1  are 
reported  the  mean  values  of  the  estimates  of  availability  at  an 
exponential  time  with  various  means  r.   Results  on  95%  confidence 
levels  are  also  given  (average  values  in  parentheses) . 

In  order  that  the  point  and  interval  estimates  behave  as  well  as 
they  do  for  small  samples  (n  ss  10)  the  bias  correction  recognizable  in 
(4.9)  is  required,  and  the  normal  percent-points  are  best  replaced  by 
Student  t  percent-points  with  n-1  degrees  of  freedom.   An  alternative 
to  the  above  procedures  is  to  apply  the  jackknife;  see  Gaver  and  Chu  . 
However,  the  present  method  is  perhaps  more  easily  carried  out  on  small 
computers  utilizing  nominal  confidence  levels  (l-a)-100%.   The  fraction 
of  the  confidence  limits  covering/surrounding  true  availability  at  r 
was  tabulated,  as  were  the  mean  upper  and  lower  confidence  limits  on 
availability.   The  following  table  provides  a  summary  of  the  results; 
as  can  be  seen  the  coverage  is  close  to  the  nominal  95%. 
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Table  5.1 

Coverage  and  Mean  Confidence  Limits,  EROBS 

Nominal  Confidence  Level  95% 

r-Val ues 

0.125  1.00  10.00 

True:     (lower)  0.92  (upper)  (lower)  0.83  (upper)  (lower)  0.80  (upper) 

Availability:  Estimated:  (0.71)  0.90  (0.96)  (0.66)  0.84  (0.93)  (0.63)  0.81  (0.91) 

(n  =  10) 
Coverage  (%)  94  96  94 

Availability:  Estimated  (0.82)  0.93  (0.97)     (0.74)  0.83  (0.90)    (0.70)  0.80  (0.88) 

(n  =  25) 
Coverage  (%)  93  97  96 


The  figures  in  parentheses  represent  mean  upper  and  lower  confidence 
limits;  notice  that,  as  anticipated,  these  tend  to  move  towards  the 
true  values  as  the  sample  size  increases. 

5.2   Predicting  Availability  under  FOBS:    Illustrations 

Again  we  use  simulation  to  assess  the  accuracy  of  the  proposed 
exponential  approximations  given  in  Section  3. 

It  has  been  found  experimentally  that  the  empi ri cal -transform- 
adapted  least-squares  approaches  of  Section  3.1  tend  to  be  numerically 
unstable  for  small  sample  sizes;  they  behave  well  when  component  models 
are  assumed  perfectly  specified  ("known")  or  for  very  large  sample 
sizes.   On  the  other  hand,  Method  3  of  Section  3.3,  adapted  to  the 
empirical  transform  as  outlined  in  Section  4.2,  performs  satisfactorily 
for  the  various  cases  considered  and  is  not  computationally  intensive. 

Simulation  results  for  three  different  estimators  for  A(t)  are 
reported.   The  estimators  are  all  of  the  form 

A(t)  =    x    +  (1  -  7r)e-/?t  (5.1) 
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where 


u  +  d 


The  estimators  differ  in  the  manner  in  which  /?  is  estimated.   The  fire 
estimate,  Ap(t)  ,  estimates  /?  by 


pp  =   (u1   +  d-1)  .  (5.; 


Note  that  Ap(t)  is  the  estimator  that  would  be  obtained  if  it  were 
simply  assumed  that  {Uj}  and  {D(}  are  independent  sequences  of 
independent  exponential  random  variables  and  maximum  likelihood  is  ugj 
to  obtain  estimates  of  A(t) ,  i.e.  if  the  simplest  Markov  model  were 
automatically  invoked.   The  second  estimator,  Ae(t),  estimates  /?  as  in 
(4.13)  with  s  =  i  .   The  third  estimator,  Aq(t)  ,  estimates  /?  as  in 
(4.15)  with  fi    =  Tjrr  .   If  the  estimated  (3 '  s  are  negative  then  /?  is  set  I 
equal  to  /?p  ,  the  MLE  estimator  of  (3     if  it  were  known  that  {Uj}  and  {Dj} 
are  independent  sequences  of  independent  identically  distributed 
exponential  random  variables.   In  our  experience  such  pathological 
cases  are  rare . 

All  simulations  were  done  on  an  IBM  3033  AP  at  the  Naval 

Postgraduate  School  using  the  LLRANDOMII  random  number  generating 

1 1 
package;  Lewis  and  Uribe   .   The  simulations  reported  in  Tables  (5.2) 


and  (5.3)  have  500  replications.   Generated  in  each  replication  is  a 
sample  of  nu  =  25  up  times  and  nd  =  25  down  times.    In  the  experiment 
reported  in  Table  (5.2)  the  up  times  and  down  times  are  independent. 
Three  distributions  are  used  to  generate  the  down  times: 
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Table  5.2 

Estimates  of  Availability 
Independent  Up  and  Down  Times 

Time      Dist     True                /?p                      Ap(t)                   0e                     Ae(t)  /?q                      Aq(t) 

A(t)         Mean     Var      Mean     MSE      Mean     Var      Mean     MSE  Mean     Var       Mean     MSE 

(Median)                                       (Median)  (Median) 

.2             A        .83             3.1         .14        .85         .03          3.8          4.2       .83         .06  3.8          4.7         .83       .06 

(3.1)                                            (3.4)  (3.4) 

.5    (Gamma)    .71             3.1         .14        .74         .05          4.6          34         .72         .06  4.5          23           .72        .06 

(3.1)                                            (3.6)  (3.7) 

1.0                       .67             3.1         .14        .68         .05           4.7          48         .68         .05  4.6          23           .68        .05 

(3.1)                                            (3.6)  (3.5) 


.2 

B 

.85 

3.1 
(3.1) 

.23 

.85 

.03 

3.4 
(3.0) 

3.2 

.85 

.05 

3.3 

(3-0) 

3.1 

.85 

.06 

.5 

(Exp) 

.74 

3.1 

(3.1) 

.23 

.74 

.04 

3.5 

(3.0) 

4.7 

.74 

.06 

3.5 

(3.0) 

4.7 

.74 

.06 

1.0 

.68 

3.1 

(3.1) 

.23 

.68 

.05 

3.7 
(3.0) 

14 

.69 

.06 

3.9 
(3-0) 

18 

.69 

.06 

.2  C      .90  4.4  6.7       .86         .06  3.4  38  .91      .03  3.5  52  .91      .03 

(3.4)  (1.7)  (1.8) 

.5     (M.  Exp.)  .87  4.4  6.7        .77  .13  2.9  57  .87      .04  2.5  14  .87      .04 

(3.4)  (1.2)  (1.2) 

1.0  .84  4.4  6.7       .72         .17  2.5  26  .82      .06  1.9  6.1  .83      .06 

(3.4)  (1.0)  (.92) 

2.0  .79  4.4  6.7        .70         .17  2.5  75  .77      .09  1.8  9.8  .78      .08 

(3.4)  (.87)  (.80) 


Table    5.3 

Estimates  oT  Availability 
Dependent  Up  and  Down  Times 

Time   Case    True         /?q             Aq(t)           /?q  Aq(t) 

A(t)      Mean    Var    Mean    MSE  Mean  Var  Mean  MSE 

0.2     DA     .94      3.1     1.7    .85     .10  1.0  .13  .94  .02 

0.5            .87      3.1     1.1     .75     .13  1.0  .07  .87  .03 

1.0            .79      3.1     1.0    .69     .11  1.0  .07  .79  .03 

2.0 .72      3.2     1.9    .67     .05  1.1  .09  -71  .02 

0.2     DB     .97      3.1     14     .91      .07  .73  5.1  .97  .02 

0.5            .93      2.4     11     .87     .07  .83  3.8  .93  .03 

1.0            .88              2.0     7.2    .83     .07  .79  .56  .88  .05 

2.0            .81      1.8     8.0    .78     .08  .83  .45  .81  .07 
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A:   P{D  >  t}  =  1  -  e"4t  -  4te"4t  ,  t  >    0  ;      (Gamma) 

B:   P{D  >  t}  =  e~2t  ,  t  >  0  ;  (Exponential) 

C:   P{D  >  t}  =  .  9e~9t  +  .le~  25t  ,  t  >    0  ;  (Mixed  Exponential) 


The  up  times  are  generated  from  an  exponential  distribution  with  unit 
mean.   The  theoretical  values  of  A(t)  in  each  of  these  cases  are 


A:   A(t)  =     .67  +  e"45t[.33  cos(1.9t)  +  .26  sin(1.9t)]   ;  (5.3) 

B:   A(t)  =  .67  +  .33e"3t  ;  (5.4) 

C:   A(t)  =  .67  +  .  24e"34t  +  .  09e"99t  .  (5.5) 


For  each  replication  of  the  simulations  whose  results  are  reported  in 
Table  (5.2)  the  three  estimates  of  A(t)  are  computed.   The  mean,  and 
mean  square  error  (MSE) ,  of  Ap(t),  Ae(t)  and  Aq(t)  are  computed:   i.e. 
the  mean  is 


-,   500  . 
Ae(t)  =  g^Q  E   Ae(k;t)  (5.6) 

k=l 


and 


1   500   -  o 

MSE  =  5(50  £   (Ae(k;t)  -  A(t))2  (5.7) 

k=l 


where  Ae(k;t)  is  the  point  estimate  at  t  in  the  kth  realization. 
Furthermore,  the  mean  and  variance  of  the  estimates  of  ft    are  displayed 

5  .  3   Discussion  of  Tables . 

Table  (5.2)  reports  results  for  simulations  with  data  sample  size 
25.  Not  surprisingly,  the  exponential  estimator  Ap(t)  has  means  equal 
to  the  values  of  A(t)  and  the  smallest  mean  square  error  in  the  case  B 
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of  exponential  down  times.   The  means  of  Ap(t)  are  somewhat  different 
from  the  A(t)'s  in  case  A  (gamma  down  times)  and  quite  different  in 
case  C  (mixed  exponential  down  times).   The  estimators  Ae(t)  and  Aq(t) 
have  means  which  are  closer  to  the  A(t)'s  than  those  of  Ap(t)  in  cases 
A  and  C.   The  means  of  Ae(t)  and  Aq(t)  are  a  bit  low  for  t  =  1,  2  in 
case  C.    In  cases  A  and  B  the  values  of  the  MSE  of  Ae(t)  and  Aq(t) 
indicate  the  greater  variability  of  these  estimators.    In  the 
exponential  case  B,  the  true  value  of  /?  is  3;  the  estimator  J3p    has 
means  closest  to  3  for  case  B  and  much  smaller  variances  compared  to  /?e 
and  /?q  .    In  the  other  two  cases  p(t)  is  not  of  the  form  e   ;  the 
theoretical  values  of  A(t)  are  given  in  (5.3)  and  (5.5) . 

In  case  A  (gamma  distributed  down  times)  the  single  parameter,  /?  , 
of  exponential  decay  is  4.5.   The  means  of  the  estimators,  /?e  and  /?q  are 
closer  to  this  value  than  those  of  /?p  .    In  case  C  of  mixed  exponential 
down  times,  the  smallest  parameter  of  exponential  decay  is  0.341.   Once 
again  the  means  of  /?e  and  /?q  are  closer  to  this  value  than  are  those  of 
/?p  .   The  variances  of  /?q  and  /?e  for  the  sample  size  of  25  can  be  large. 
In  many  cases  the  variance  of  /?q  is  less  than  that  for  (3e  .    Increasing 
the  sample  size  decreases  the  variances,  as  is  anticipated. 

Table  (5.3)  summarizes  a  simulation  study  of  the  procedure  of 
Section  4.2  in  two  cases  in  which  the  pairs  (Uj,Dj)  are  independently 
and  identically  distributed  with  Uj  and  Dj  perfectly  dependent.   The  two 
cases  are : 

DA:   Uj  =  Ej   and   Dj  =  ^Ej   with  {Ej}  independent  identically 
distributed  exponentials  with  unit  mean. 


DB:   Uj  =  Ej  and  Dj  = 


'^Ej   with  probability  .9, 
,4E:   with  probability  .1 
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with  {Ej}  independent:  exponential  random  variables  with  unit  mean.    In 
case  DA 


A(t)  = 


2  +  I  e"4 

3  +  3 


(5.8 


In  case  DB 


A(t)  =  2  +  Q.lle"1  +  0.22e-°-27t 


(5 


In  the  simulations  two  estimators  for  Fc(s)  are  used.   The  first 
(incorrectly)  assumes  that  {Uj}-  and  {Dj}  are  independent  and  uses 
<f)y  (s)  <f>D  (s)  to  estimate  Fc(s)  .   The  second  does  not  assume  independence, 
and  instead  estimates  Fc(s)  by 


*c(s)   =  i  E 


!   5c   -stuj+dj] 


(5.10 


i=l 


The  table  shows  results  for  the  estimators  Aq(t)  and  /?q  .   The  sample 
size  is  25  and  there  are  500  replications. 

Table  (5.3)  informs  us  that  the  estimators  may  be  noticeably 
sensitive  to  the  choice  of  the  estimator  of  Fc(s) .   If  Uj  and  Dj  are 
dependent,  then  using  ^y(s)</>D(s)  to  estimate  Fc(s)  can  be  quite 
misleading.   As  a  result,  it  is  suggested  that  if  there  is  a 
possibility  that  Uj  and  D:  may  be  dependent,  ^c(s)  is  the  more  model- 
robust  estimate  of  Fc(s),  and  hence  of  the  desired  availabilities. 

Note  that  in  case  DA,  A(t)  has  exactly  the  presumed  form  of 
(4.10) .   When  0C  is  used  in  this  case,  the  mean  value  of  the  0q ' s  equal 
the  theoretical  value  of  unity  in  all  but  the  case  t  =  2.0.   In  case 
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DB ,  p(t)  is  actually  a  mixture  of  two  exponentials,  one  having  rate 
unity  and  the  other  having  rate  0.27;  when  <f>Q     is  being  used,  the  mean 
values  of  the  /3q  '  s  fall  between  the  two  correct  rates.   Emphasis  of 
either  rate,  and  the  corresponding  probability,  can  be  achieved  by 
we  ight  i  ng . 

5 . 4   Semi -Parametr i c  Confidence  Intervals  Illustrated . 

Table  (5.4)  shows  bootstrap  confidence  intervals  for  Aq(t)  and  /?q 
at  various  times.   A  single  sample  of  25  up  times  and  25  down  times  is 
generated  by  simulation.   The  up  times  are  independent  with  unit-mean 
exponential  distribution;  they  are  independent  of  the  down  times.   The 
down  times  are  independent  with  the  gamma  distribution,  case  B.   One 
hundred  bootstrap  replications  were  then  carried  out  and  Aq  and  /?q  are 
computed  for  each  replication.   Both  estimators  of  Fc(s)  ,  ^y  (s)  <£D  (s) 
and  <^c(s),  are  used.   Table  (5.4)  exhibits  the  5th,  50th,  and  95th 
order  statistics  of  the  estimates  which  give  95%  confidence  intervals 
of  the  parameters. 

All  confidence  intervals  for  (3q    cover  the  true  values  for  the 
particular  sample  utilized.   The  confidence  intervals  for  f3q    cover  the 
exponential  decay  parameter  of  value  4.5,  for  the  gamma  case.   The 
confidence  intervals  for  Aq(t)  are  the  same  for  both  methods  of 
estimating  Fc(s)  suggesting  that  not  much  is  lost  by  using  </>c(s)  to 
estimate  Fc(s)  in  the  independent  case.   As  a  result  it  is  suggested 
that  if  there  is  a  possibility  that  {Uj  and  {D(}    may  be  dependent  then 
<£c(s)  is  the  more  robust  estimate  of  Fc(s). 

A  similar  bootstrap  experiment  was  carried  out  with  a  single 
sample  of  25  up  times  and  25  down  times  generated  from  model  DA  in 
which  the  up  and  down  times  are  dependent.   One  hundred  bootstrap 
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Table  5.4 

Bootstrap  Confidence  Intervals 

Gamma  Distributed  Down  Times 

Independent  Up  and  Down  Times 


Time          True                 0q  Aq(t)  f3q  Aq(t) 

A(t)  .05  .50     .95  .05  .50     .95  .05  .50     .95  .05  .50     .95 

.2              .83  1.6  4.1    13.6  .71  .84     .92  1.6  3.9    11.8  .71  .84     .92 

.5              .71  2.4  4.1    12.9  .61  .73     .82  2.4  3.9     9.8  .61  .73     .82 

1              .67  2.3  3.4    16.3  .58  .71     .79  2.2  3.4    19.1  .58  .71     .79 


replications  were  generated  and  Aq  and  /?q  computed  for  each  replication 
with  both  estimators  of  Fc(s).   Once  again  the  resulting  confidence 
intervals  for  A(t)  were  very  similar  for  both  estimators  of  F<-(s)  . 

Of  course  the  above  results  are  quite  fragmentary,  but  seem  useful 
and  promising.   Further  sampling  experiments  and  asymptotic  analyses 
will  shed  more  light  on  the  behavior  of  the  estimating  procedures 
explored,  and  may  well  suggest  alternations  or  replacements. 
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6.   CONCLUDING  COMMENTS. 

This  paper  argues  that  probability  models  of  the  availability  of 
various  systems  can  be  expressed  in  terms  of  Laplace  transforms,  and 
that  the  finite-time  behavior  of  such  systems  can  be  inferred,  nearly 
non-parametr i cal ly ,  from  data.   Our  approach  has  been  to  invoke  the 
empirical  Laplace  transforms  and  to  utilize  its  easy  direct  interpreta- 
tion (EROBS) ,  in  conjunction  with  a  presumed  approximate  exponential 
rate  of  approach  to  steady  state,  to  deduce  availability  at  a  fixed 
finite  time  after  a  known  initial  moment  (FOBS). 

The  methods  proposed  are  distinguished  by  their  simplicity  and 
moderate  computer  intensivity  as  well  as  by  their  lack  of  direct 
dependence  upon  probability  models  in  "up"  and  "down"  times  selected 
from  conventional  families  such  as  the  Gamma.   Simulations  have  been 
used  to  evaluate  the  procedures  suggested,  and  to  provide  approximate 
confidence  limits,  either  by  asumptotics  (utilizing  the  approximately 
Normal /Gauss i an  behavior  of  the  empirical  transform),  or  by  a  simple 
re-sampling,  Efron's  bootstrap   .    In  particular,  we  have  examined  the 
effect  on  inference  quality  of  assuming  the  wrong  joint  probability 
model:   one  that  falaciously  assumes  independence  when  dependence 
(between  up  and  down  times)  is  actually  present. 

No  claim  is  made  that  the  methods  proposed  here  are  the  best 
available;  in  fact  there  are  many  alternatives.   One  is  to  analytically 
invert  the  empirical  transform  of  availability  possibly  by  use  of  the 
Stehfert  algorithm   ,  although  competitors  are  available.   Another  is 
to  bootstrap  directly;  the  latter  exercise  involves  re-sampling  up  and 
down  times  from  the  observed  data,  to  reconstruct  the  sample  path  of 
the  process,  and  to  score  1  at  time  t  if  an  up  time  covers  t,  otherwise 
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score  0.   Finally,  A(t)  is  estimated  by  "the  proportion  of  (re)samples 
that  count  1.   Confidence  limits  are  available  from  the  basic  bootstra 
technology.   The  direct  bootstrap  approach  is  being  investigated  by 
Lee    .   Another  option  is  to  replace  the  distribution  functions  in  the 
renewal  equations  for  availability  by  their  empirical  counterparts  and 
numerically  solve  the  empirical  renewal  (Volterra-type)  equations,  wit 
subsequent  bootstrap  follow-ups  to  assess  uncertainty.   All  such 
methods  promise  to  be  far  more  computationally  intensive  than  our 
present  approximate  approaches.   Their  investigation  has  been  deferred 
Application  of  our  approximation  procedure  to  infer  the  M/G/l  queue 


finite-time  behavior  is  under  way;  Jacobs  and  Gaver 


11 
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